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ABSTRACT 

We use Schwarzschild's orbit-superposition technique to construct self-consistent mod- 
els of the Galactic bar. Using minimisation, we find that the best-fit Galactic bar 
model has a pattern speed fip = 60 km s~^ kpc^^, disk mass Md = 1.0 x IO^^Mq 
and bar angle ^bar — 20° for an adopted bar mass Mbar — "2 x 1O^°M0. The model 
can reproduce not only the three-dimensional and projected density distributions but 
also velocity and velocity dispersion data from the BRAVA survey. We also predict 
the proper motions in the range I = [—12°, 12°], b = [—10°, 10°], which appear to be 
higher than observations in the longitudinal direction. The model is stable within a 
timescale of 0.5 Gyr, but appears to deviate from steady-state on longer timescales. 
Our model can be further tested by future observations such as those from GAIA. 

Key words: Galaxy: bulge - Galaxy: centre - Galaxy: structure - Galaxy: general - 
Galaxy: nucleus - Galaxy: formation 



1 INTRODUCTION 

It is well known that most s piral galcixies hos t a bar structure 
in their central region (e.g. lLee etal] 120111 ). Therefore, one 
of the most important issues in galaxy formation and evo- 
lution is to understand the structure and dynamical prop- 
erties of barred galaxies. Our own Milky Way is the near- 
est barred galaxy. Compared with other distant galaxies, 
the Galaxy has extensive observed photometric and kine- 
matic data which enable us to study the bar structure in 
detail. A thorough understanding of the structure and dy- 
namics of the Galactic bar may help us understand the for- 
mation of other spiral galaxies, and test the validity of the 
popular Cold Da rk Matter structure formation model (e.g. 
IShen et al.|[201ol ). 

Observationally, the Galactic bar model can be con- 
strained by surface b rightn e ss maps jDwek et al.l 1 19951 : 
iBabusiaux fc Gilmord l2005l : IConzal ez et al.' '201 ih. mi- 
crolensing optica l depth maps ( Blitz fc Spergel 



Zhao et al.l 19951) and star counts (jStanek et al.l 



1991 



1994 



Mao fc Paczvhskil |2002| : iRattenburv et al.1 120071 ). Different 



observational techniques, wavelengths and fields probe dif- 
ferent aspe cts of the bar. Schwar zschild's orbit-superposition 
technique (jSchwarzschildl [l979l ) provides us the possibility 
to construct a model that can fit all observations. Using the 



orbit-superposition technique, Zhao (1996, hereafter ZH96) 
constructed a self-consistent Galactic bar model, which can 
fit the surface brightness, velocity and velocity dispersion 
in the Baade's window (BW). However, recently it was 
found that the predicted rotation curve in ZH96 is in- 
consistent with the results from the Bulge Radial Veloc- 
ity Assay (BRAVA , iRich et al.l l2007l : [Howard et all l2008l : 
iKunder et al]|2012l) . Possible explanations are (1) the data 
resolution in ZH96 is not sufficiently high, and/or (2) the 
initial conditions of orbits include too many loop orbits. 
Also using the orbit-superposition technique, iHafner et al.l 
(2000) constructed a dynamical model of the inner Galaxy, 
which fits most of the available data in one very inten- 
sively observed bulge field, Baade's Window {I, b) = (0, —4); 
the only disadvantage is that the proper motion disper- 
sion in their mod el is higher than the measured values by 
ISpaenhauer et aL 

I [l99Z) . Another useful method to gener- 
ate self-con sistent dynamical model s is the made-to-measure 
algorithm (ISver fc Trema inj 'l996|; Ide Lorenzi et al.l l2007l : 
lDehnen|[2009l : iLong fc^ ao 2010), which was implemented 
for the Milk Way by Bissan tz et afl (|2004l ). However, only 
the density map of Milk Way was used to construct the dy- 
namical model, no kinematic constraints were applied and 
their effective field is small. Furthermore, the effective par- 
ticle num ber in their study turns o ut to be small (only a few 
thousand. IRattenburv et al]|2007l ). 
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Recently, extensive observations of the central re- 
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line of sight 
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Figure 1. The coordinate system used throughout the paper. 
The major axis (x-axis) of the bar is misaligned with the line of 
sight (X-axis) with an angle 9bar- The z-axis comes out of the 
page. The bar rotation and the rotation of the Sun with respect 
to the Milky Way centre are indicated. Notice the rotation of 
the bar is in the negative z-axis and the near-side of the bar has 
positive longitude (1). 



gion of the Milky W a y by 
scope (iKozlowsk i et al.] 
scopes (BRA VA: iRich et al 



2006h 



the Hubble Space Tele- 
an d ground-based tele- 
20071: iHoward et all |200S" 



Kunder et ai]|2012l : OGLE: lUdalski et al.ll2000l : ISumi et all 
2004 ) provide many other large samples of kinematic data. 
In this paper, our aim is to construct a self-consistent and 
stable bar model which can fit all the presently-available ob- 
served data in the central region of the Milky Way by using 
the Schwarzschild method. As we will see later, we are able 
to produce a self-consistent model, but there are issues with 
long term stability. 

The paper is organized as follows. In section 2, we de- 
scribe the density and potential model for this study. Sec- 
tion 3 presents the details of our implementation of the 
Schwarzschild method. In section 4, we show the main re- 
sults. Section 5 shows the stability of the bar model. Con- 
clusion and discussion are given in section 6. Throughout 
this paper, we adopt the distance to the Galactic center as 
J?o = 8 kpc. 



2 DENSITY AND POTENTIAL OF THE 
CENTRAL REGION OF THE GALAXY 

For clarity, in Figure [U we first establish the coordinate sys- 
tem we use throughout this paper. In particular, the major 
axis of the bar is along the x-axis, and is at an angle ^bar 
with respective to the line of sight (the X-axis). Notice that 
the bar angular momentum is along the negative z-axis, but 
we still write the pattern speed as positive for abbreviation. 



used the star counts of the red clump giant stars to con- 
strain the triaxial Galactic bar mo dels. Both IStanek et al.l 
l|l997h and iRattenburv et all ||2007^ found the analytic bar 
model given in iDwek et al.ni995l) fit the data well. In this 
paper, we also adopt the Dwek et al.l (|l995l ) bar model as in 
ZH96, which has the following form 



p{x,y,z) = po 



exp 



''l 1 I -1.85 



exp(- 



r2 



(1) 



where the first term represents a bar with a Gaussian ra- 
dial profile and the second term a spheroidal nucleus with a 
steep inner power law and an exponential outer profile. In 
this paper, we refer both the bar and nuclear components 
simply as the "bar" . The central density po is determined by 
normalising the total mass of the bar, Mbar = 2.0x 10^°Mq, 
which is fixed throughout the paper. The radial functions ri 
and r2 are defined as 



and 



r2 



+ 



+ 



1/4 



q^{x^ + y'^) + z'^ 



1/2 



(2) 



(3) 



where the principal axes of the bar are xq — 1.49 kpc, j/o ~ 
0.58 kpc, zo = 0.40 kpc and the bulge axis ratio is q — 0.6. 



2.2 Disc potential and density 

We do not include explicitly any dark halo in our potential. 
iKlypin et all l|2002l ) shows that the dark halo has to be very 
low in mass in the central part in order to allow many mi- 
crolensing events by baryonic material. Instead of the usual 
halo plus exponential disk, we use a Miyamoto-Nagai po- 
tential to represent the disk plus halo. The disk potential is 
given by 



'^d{x,y,z) 



where 



2 , 2 , 

X +y + 



IMN + {z + &mn) 



1/2 



(4) 



(5) 



flMN = 6.5 kpc, 6mn = 0.26 kpc and Md is the total disc 
mass. In ZH06, Afd = SAfbar and Afbar = 2 x 10^° A/© is the 
total mass of the bar. In this paper, we will consider different 
values of the disk mass but keep the bar mass fixed. 

The density distribution of the MN disk is given by 



2.1 Density and potential of the bar and bulge 

Lots of observation data have been used to constrain dy- 
namical models of the Galaxy. Using the COBE D IRBE 
multiwavelength observations, iDwek et al.l |[r995l) con - 
structed many analyt i c bar mod els. IStanek et al.l lll997h . 
iBabusiaux fc Gilmord l|2005h and IRattenburv et al.l (|2007l ') 



Pd{x,y,z) = 



flMN-R^ + (flMN + 3v^2M~62^)(aMN + z"^ + &mn)^ 
[R^ + {aUN + ^Z^ + hi,^Y?/\^2 + f,2^^)3/2 

2 I 2 

^x +y . 
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Table 1. Twenty Hernquist-Ostriker expansion coefficients A„ 
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Figure 2. Density distribution along the major (a;-) axis. The 
solid, dotted and dashed lines represent the bar, disk and total 
density of the model, respectively. 



2.3 Density, potential and accelerations of the 
system 

Figure [2] shows the density distribution for the models along 
the major axis (a;-axis). It is obvious that the bar dominates 
the mass distribution of the system in the inner 3 kpc. 

In order to solve the Pois son's equat ion, we follow 
iHernauist fc OstrikeJ (|l992l ) and lZhad (|l996l ) to expand the 
potential and density on a set of simple orthogonal basis of 
potential-density pairs in the spherical coordinates 

4>n!m, (6) 

^= I 

n,i ,m 



"l>bar = 

where 



(1 + S)2' + 1 

r-s = 1 kpc 



g: 



(21+3/2) I S 



Pim (cos 61) cos(m<ji), (7) 
(8) 



and G^'^^'''^^^^) is the Gegenbauer polynomial of ^. The 
expansion coefBcient A„im can be calculated by 



RR{r)Pim{X)cos{m.4>) 



1 r2TT 

dX / d(l>p{r,e,(j)) 
1 Jo 



(21+3/2) 



(9) 



(10) 



where £ = ^-r^ , ^ = cos 

r _ r. ^ r(n + 4/ + 3) 

- ^"'2«+6n!(n + 2/ + 3/2)[r(2/ + 3/2)]2^ + """^ 



2 (/ + m)! 

21 + 1 {I -my: 



(11) 



K^i = -n{n + + i) + [l + 1){21 + 1) (12) 

and Srno is the Dirac Delta function which is defined as Smo = 
1 for m = Q and SmO = for m 7^ 0. 
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Each expansion coefficient is determined by the three 
quantum numbers n, I, and m. For a triaxial model, only 
the even quantum number terms are non-zero. The circular 
velocity profile is shown in Figure [3] along the intermedi- 
ate axis (y-axis). As can be seen, there is little difference 
between models with 20 and 40 terms (left panel). So in 
this paper we will adopt twenty expansion coefficients in 
the orbit integration to save CPU time. They are listed in 
Table [T] The right panel of Figure shows the dependence 
of the circular velocity on the disk mass. Clearly, the model 
with the massive disk mass has the large circular velocity 
beyond 1 kpc. Figure 0] shows the contour of effective po- 
tentials of the model in the x — y plane (top panel) and x — z 
planes (bottom panel). 

Since we have the potential for the systems, the accel- 
erations can be easily calculated from the potential; we do 
not give detailed expressions for the acceleration here. 



3 MODEL CONSTRUCTION 

3.1 Orbit-superposition technology 

bmce ISchwarzschiidI (| 19791 ') pioneered the orbit- 
superposition method to construct self-consistent models 
for three-dimensional mass distribution, it has been widely 



applied in dyn amical modelling 
van der Marel et al. 1998: Binncy 



(e.g. 

2005; 



Rix et al. 1997: 



van de Ven et al 



Capuzzo-Dolcetta et al.l l2007l: Ivan den Bosch et al 



200a, , - ■ ■ r 

12008 : I Wang et all I2OO8I : IWu et al.l 120091 1. The kev point 
of this method is to construct an orbit library which 
is sufficiently comprehensive in order to reproduce the 
available observations. 

Specifically, let A^o be the total number of orbits and 
A'c be the number of spatial cells. For each orbit j, we count 
the fraction time Oij and projected quantities Pij that they 
spend in each cell i. The fraction time of Oij of each orbit is 
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Figure 3. Circular velocity along the intermediate (y-) axis of 
the bar. Left: The solid, dotted and dashed lines represent the 
results from 10, 20 and 40 terms of the Hernquist-Ostriker expan- 
sions, respectively. The disk mass is Md = 1.0 X 1O"M0. The flat 
horizontal line indicates an amplitude of 200kms~^. Right: The 
solid, dotted and dashed lines represent results of = 2.25, 
5 and SMb^r, respectively. Twenty terms of Hernquist-Ostriker 
expansions are used. 



obtained as follows. Every orbit is integrated for one Hubble 
time (tfr), and TV — 10,000 output (position and velocity) 
are stored at constant time interval for each each orbit. If 
the orbit j crosses the cell i one time, we will increase the 
number Nij by 1. Then, the fraction Oij is determined by 
Oij = Nij/N. The orbit weight Wj for each orbit j is then 
determined by the following equations: 



(13) 



where jii can be the volume density, surface density, or mo- 
ments of the velocity distribution in each cell i. If is the 
meiss or the veloc ity, then equation l|13p is same as equation 
(2.4) or (2.5) of iPfennigerl (|l984 ). Following ZH96, we di- 
vide the first octant into 1000 cells with similar masses. Due 
to symmetry of the model, the other octants are "reflected" 
to the first octant. In each (x-, y-, and z-) direction, the 
system is divided into 10 bins. Each cell is a small box with 
dx = 0.25 kpc, dy = 0.15 kpc and dz = 0.10 kpc. The central 
cell covers the box with x = [-0.25, 0.25], y = [-0.15, 0.15] 
and 2 = [-0.1,0.1]. 

More practically, equation [T3] can be written as a set of 
linear equations 



(14) 



We adopt the no n-negative least squares (NNLS) 
method (IPfennigeij 1 19841 ) to solve Wj, i.e., the following Xmj 



1=1 



(15) 



is minimized with respect to Wj {j = 1, , A'^o) to obtain 

the values of Wj ^ 0. It is obvious that the NNLS fit will 
find a unique solution if the number of orbits is smaller than 
the number of constraints. However, such a model may not 
be self-consistent. A more meaningful result with the NNLS 



u 

CL 




u 

CL 




Figure 4. Top: Effective-potential contours for 20 shells with 
equal mass in the x—y plane. Bottom: Effective-potential contours 
in the x — z plane. The value of effective potentials in the most in- 
ner and outer surfaces are —1.81 X 10^ and — l.OSx 10^(kms~^)'^, 
respectively. 



method should use a large number of smooth orbits well 
sampled in the phase space. In this case, many exact solu- 
tions with = are possible. The NNLS method will select 
one of the possible solutions. 

Our density model is smooth, therefore, we expect 
that the phase-space density from the reconstructed self- 
consistent model is also smooth. We employ a simple 
smoothing procedure to fit the data. We require simply that 
the orbits wi th adjacent initial condit ions have nearly the 
same weig ht l|Merritt fc Fridman|[T99^ . hereafter MF96). In 
this approach, Equation [T5l becomes. 



2 _ 



°ij ) Oij V V J 



(16) 



j=i 



where A is a smoothing parameter and A = A'^^. 

MF96's method strongly depends on the size of the cells. 
In order to assess whether the smoothing method affects 
our results, we also used another smoothing method, also 
adopted by ZH96. The key point of this smoothing method is 
that orbits with similar integrals of motion should have sim- 
ilar orbit weights. In a rotational bar system, only Jacobi's 
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Figure 5. Radial velocity and velocity dispersion along four strips 
with b= -4°,b = -6°, b = -8° and 1 = 0°. The solid (SMI) and 
dotted (SM2) lines indicate the results from our models using the 
smoothing method as shown in M96 and ZH96, res pectively. The 
filled circles with error bars are data from BRAVA ijKunder et alj 
I2OI2I) . 



integ ral Ej is an integral of motion ijBinnev fc Tremaind 
Il987h . here we take the time averaged quantities (Lz) and 
(L^) as the effective integrals, where and are an or- 
bit's instantaneous angular momentum components along 
the minor and major axes, respectively. Orbits with simi- 
lar Ej, (Lz) and (L^) have similar orbit weights. We have 
checked the results from these two smoothing methods and 
found no significant difference (See Figure [5]) in their abili- 
ties to fit the BRAVA data. The only difference is that we 
find the MF96's method yields more orbits with non-zero 
weights than those of ZH96. From now on, we only present 
results using the smoothing method in M96. 



3.2 Constructing the orbit library 

In a rotational triaxial system, as mentioned before only 
Jacobi's energy is an integral of moti on. Early studies 
have shown that most orbits are chaotic (|VogUs et al.ll2007l : 
iManos fc Athanassoulalboill . see also sectionl4.3.2p and the 
typical regular orbits are xi -type (|Binnev fc Tremainell 19871 : 
IZhaoll994l : IZhao et al]|l 994'') in a bar model. We do not know 
the explicit phase space distribution f{x,y,z,Vx,Vy,Vz) for 
the chaotic orbits because they lack integrals of motion in 
a rotating bar potential, we will thus consider two differ- 
ent methods to generate the initial conditions of orbits. The 
common point of these methods is that the bar model is di- 
vided into 20 shells with nearly equal mass spatially along 
the a;-axis by using the Monte Carlo integration (see Fig. Q. 

The first method (ICl) is similar to the one adopted in 
ZH96. Here we give a brief description of its main ingredi- 
ents, and refer the reader to Appendix B in ZH96 for more 
details. The orbits are launched in close pairs perpendicular 
to the xz~, yz- or xy— symmetry plane or the x, y axis. 
The initial position of each orbit in each shell has the same 
effective potential, which is defined as that in the bounded 
surface of two close shells. The initial velocity is tangential 
and less than the local circular velocity. Only 1000 initial 
conditions of orbits were generated in ZH96, here we extend 
the number of orbits to ~ 20, 000. 
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Figure 6. Density contours for the input model (solid lines) and 
from orbits (dashed lines) in the x—y and y—z planes. The left and 
right panels represent the results from ICl and IC2, respectively 
(see 



The second meth od (IC2) is similar to the one used 
in lHafner et al.l (|2000 '). The initial conditions are generated 
with known distribution functions / = X]i=i Cip{x,y, z)hi, 
where Ci,i = 1,2,3 are three normalizing constants, p is 
the density distribution of the system. The three functions 
hi-1,2,3 are defined as 



h3{Vx,Vy,Vz 

and 



(2^)3/V. 



■ exp 



2(72 2a| 



(17) 



hl,2{vR, V^,Vz) = 



1 



X exp 



A. 



2gI 



2ol 



where Vd is defined by 



Vci{R) = 250 km s" 



1 + 



0.1 kpc 
R 



(18) 



(19) 



For the velocity dispersion par ameters, w e adop t the same 
values as listed in Table 2 of H afner et"al] l|2000l ). We select 
50,000 initial conditions using this method. 

In Figures |S] and [T] we present the volume density and 
the projected density contours for model 25 (see Table 1), 
respectively. In each figure, the solid and dashed lines rep- 
resent the results from the input and orbit models, respec- 
tively. The left and right panels represent the results from 
the two methods (ICl and IC2), respectively. It is seen that 
the orbit from both methods can reproduce the density dis- 
tribution and projected density distribution well. The differ- 
ence between the reconstructed and input densities is small, 
which indicates that our model is self-consistent. 

Figure [8] further compares the velocity and velocity dis- 
persion from the model with those from the BRAVA data. 
There is no significant difference between the results from 
these two initial conditions, thus from now on, we only 
present results from th e first method (I Cl). 

As pointed out bv lPfenniged (|l984l ). a reasonable orbit 
integration time may be determined by the fluctuation of the 
Oij between two successive halves of the integration time. 
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Figure 7. Similar to Figure (6] but for the projected surface den- 
sity. 
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Figure 8. BRAVA data versus model as in Figure [5] The solid 
(ICl) and dashed (IC2) lines indicate the results from our model, 
while the filled circles with error bars are data from BRAVA. 
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Figure 9. Mass distribution for the input model (solid lines) and 
from orbits (dashed lines) in the x—y and y—z planes. The left and 
right panels represent the results from the first and second half 
Hubble time. Top: Results for using 1000 cells. Bottom: Results 
for using 400 cells. 



The Oij actually reflects the orbit densities and a superpo- 
sition of Oij reflects the system density. For regular orbits, 
Oij 's can reach stable values in a relatively short time. How- 
ever, Oij 's for irregular orbits only can converge after a very 
long integration tim e, at least 1000 Hubble times as sug- 
gested bv Pfennigeil (|l984l ). In practice, it will be too time- 
consuming to integrate a large number of orbits for such a 
long time. We have compared the fluctuation of the Oij 's be- 
tween two successive halves in one Hubble time with those 
in ten Hubble times for irregular orbits, there is no clear 
improvement in the convergence. In this paper, the orbits 
are integrated for one Hubble time unless stated otherwise. 
Although the typical fluctuation of Oij 's is about ~ 20% be- 
tween the flrst and second half, the typical mass fluctuation 
in each cell is small (below 2%) , as can be seen from Figured 
Figure [9] shows the mass distribution when the integration 
time is the first and second half Hubble time respectively. 
We only consider the mass distribution in Equation psp . 
and no smoothing is adopted. In order to decrease the fluc- 
tuation, we decrease the spatial resolution of cells, reducing 
the cell number from 1000 to 400 by merging every two and 
half adjacent cells into one along the z-direction. The typi- 
cal fluctuation of the On between the first and second half 



Hubble time is reduced to ~ 10% and the mass fluctuation 
is ~ 0.05%. 

Since there are many parameters in the bar model, un- 
less stated otherwise, we adopt Mbar = 2.0 X 10^° Mq, Md = 
SMbar, f2p = 60 km s~^ kpc~^ and bar angle ^bar ~ 13.4°. 
These parameters are the same as those used in ZH96 (model 
25 in Table H} for ease of comparisons with ZH96. 



4 RESULTS 

4.1 Model constraints 

The key point is to solve the orbit weights from equation 
1151 To do this, the volume density, projected density, radial 
velocity and velocity dispersion along four windows b = —4°, 
b = —6°, b — —8° and I — 0° (see Figure (5)1 are used as 
constraints to solve the weight of each orbit. The volume 
density and projected density are obtained directly from the 
density distribution of the bar and Miyamoto-Nagai disk 
model. Since our aim is to construct a self-consistent bar 
model, the volume density is fitted only inner 3kpc around 
the Galactic center. 
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Figure 10. Observed fields by the BRAVA survey from 2005 to 
2008 (also see Fig. 1 in Kunder et al. (2012)). 



The kinematic constraints (radial velocity and disper- 
sion) are from the BRAVA survey, conducted from 2005 to 
2008 in 8746 fields, which are shown in Figure [TOl It is seen 
that most data are located in the range / = [—12°, 12°] and 
b = [—10°, 10°], where I is the Galactic longitude and b is 
the Galactic latitude. In this paper, we fit the projected den- 
sity, the radial velocity and velocity dispersion in the range 
1 = [-12°, 12°], b= [-10°, 10°]. 

4.2 Dependence on model parameters 

In this subsection, we vary the model parameters around the 
ZII96 model (model 25 in Table[2l) to see the trends. We will 
explore the parameter space more systematically in section 
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Figure 11. Projected density of the input bar model for differ- 
ent bar angles. The solid, dotted, and dashed lines represent the 
results for the bar angle 13.4°, 20°, and 30°, respectively. 
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4. 2.1 Bar angle 

From the COBE map, we know that there is a clear offset 
between the major axis of the bar and the line-of-sight to the 
Galactic center. However, it s value is not accura tely known. 
From COBE observations, iDwek et all (|l995l ') found the 
value of bar angle is 20° ± 10°. Using the sam e COBE map, 
IZhad I 1994f ) obtai ned a bar ang l e 13.4 ° while [sinnev et al] 
l|l997t) found 20°. lAlcock et alT ('2OOO') found this value to 
be 15°. Recently, the 6.7-G Hz methanol mas ers showed a 
45° orientation of bar ang le (|Green et al.ll201ll ). No consen- 
sus appears to be emerging among the recent observations 
concerning the bar angle. 

In Figure [TT] we show the projected density maps of the 
input bar model for different bar angle. It is clear that the 
isodensity maps become more sharp-edged with an increas- 
ing bar angle, which means that the bar angle can be deter- 
mined if high quality surface brightness data are available. 
In this paper, our bar model is from ZH96, which attempts 
to match COBE observations. We can vary the bar angle 
from 13.4° to 20° or even 30° within the error bar of the 
model. However, a large bar angle (40°) would require us to 
refit data in terms of other bar parameters (lengths and ax- 
ial ratios. [Zhao et al]|l996l V Therefore, we restrict ourselves 
to three different bar angles 13.4°, 20° and 30°. 

Figure [12] compares the radial velocity and velocity dis- 
persion from the orbit projection with ones from the BRAVA 



Figure 12. Radial velocity and velocity dispersion along the 
bulge major axis for different bar angles. The pattern speed 
is f2p = 60 kms~^kpc~^ and = SMban where Mb^r = 

2 X 101°Mq. 

data for different angles. There is a small difference between 
the results from different bar angles. However, it is clear that 
kinematics alone constrain the bar angle poorly. 

4-2.2 Pattern speed 

The pattern speed of the Galactic bar has been esti- 
mated from diffe r ent m ethods and are somewhat uncertain. 
iDebattista et al] (|2002l ) used the Tremaine- Weinberg con- 
tinuity method to the OH/IR stars and obtained a value 

= (59 ± 5 ± 10)kms-^ kpc'^ lEnglmaier fc Gerhard! 
(|l999t ) obtained 57p « 60kms ^kpc ^ by comparing 
the gas flow in hydrodynamic simulations with the ve- 
locity curve fro m HI and CO observations. From the 
length of the bar llBinnev et aL l ll997l : iBeniamin et aLll2005l : 
ICaijrera-Lavers et al.ll2007l ). the pattern speed was given in 
a wide range fip ~ (35 - 60) kms^^kpc"^ (Gerhard 2010) . 

We consider four pattern speeds fip = 40, 50, 60 and 
80 kms~^kpc~^. In Figure [Ol we present the velocity and 
velocity dispersion distributions for different pattern speeds. 
Obviously, the velocity dispersion profile strongly depends 
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Figure 13. Similar to Figure [T2l but for different pattern speeds 
of the bar. The soUd, dotted dashed and long dashed lines rep- 
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Figure 14. Similar to Figure [T2l but for different values of the 
disk mass. The solid, dotted and dashed curves represent = 
2.25, 5 and 8 Afban respectively. f!p = 60 km s~^ kpc~^, S^ar = 
20°. 



on the pattern speed. The predicted velocity dispersion of a 
model is inversely correlated with its pattern speed. 



4-2.3 Disc mass 

The disk mass is another parameter which is not accurately 
known. We consider three different values of the disk mass. 
Figure [14] shows the dependence of the velocity and velocity 
dispersion on the disk mass. As can been seen, the veloc- 
ity dispersion profile from the model strongly depends on 
the value of the disk mass: as expected, a less massive disk 
induces a lower velocity dispersion than a more massive one. 



4.3 Best-fit model 

As mentioned in section \A72\ the kinematics from the model 
depend on the bar angle, pattern speed and disk mass. In 
principle, we can divide the parameter space of the bar angle, 
patter speed and disk mass into many cells. In each param- 
eter cell, we can run the orbit-superposition technique and 



use the fit to find the best-fit parameters. However, nu- 
merical calculation is expensive. Here, we calculate 36 mod- 
els with different parameters (bar angle, pattern speed and 
disk mass). The is defined as 



X 



(20) 



where Nobs is the total number of observed data, yahs 
and 1/modoi are the observed and model kinematics, respec- 
tively. Since the three-dimensional and projected densities 
are given by the input model, we do not include them in the 
fitting, although we do compare the predicted distribu- 
tions with data by eyes. 

Table [5] lists the value of for fitting velocity and 
velocity dispersion along both the major and minor axis of 
the bar for 36 models. It is seen that models 12, 13, 14, 
15, 22, 23, 24, 33 have smaller values of than others, 
visual examination indicates that model 23 is the best-fit 
model in both the three-dimensional and projected density 
distributions. Therefore, we choose it as the best-fit model 
for the BRAVA data. 

In Figures[T5]and 1161 we present the volume density and 
the projected density contours, respectively. In each figure, 
the solid and dashed lines represent the results from the 
input and orbit models, respectively. It is seen that the orbit 
from model 23 can reproduce the density distribution and 
projected density distribution well. Moreover, we also use a 
parameter 5 to describe the departure from self-consistency 
for model 23, which is defined as (MF96): 



5 = n/^/M, 



(21) 



where M is the average mass in each cell, if the total mass 
is 1, then M = 1/Nc- The is obtained from equation [15] 
by only using the mass constraints and without smoothing. 
Figure [17] shows the departure from self-consistency as a 
function of the number of orbits. It is noted that departure 
parameter 5 strongly depends on the number of orbits and 
the cell number, 5 is smaller than 10~® (0.013) when 17,323 
orbits are adopted in equation 1151 by using the 400 (1000) 
mass cells, which again shows that model 23 is nearly self- 
consistent. 

The velocity and velocity dispersion distributions of 
model 23 in Figure 1181 Note that model 23 can fit the 
BRAVA data well with only a few outliers. 



4-3.1 Predicted proper motions 

Proper motions are not taken into consideration in solving 
Equations [15] and [16] because the absolute values of proper 
motions can not be obtained from observations at present 
due to the lack of absolute astrometry. We can still com- 
pare the predicted proper motion dispersions with those 
observed. Table |3] shows the proper motions in some ob- 
served fields together with the predictions from model 23. 
It is seen that the proper motions along the latitude from 
model 23 in the Baade's and Sagittarius's Window are in 
good agreement with those in observations, while the lati- 
tudinal proper motion in the Plaut's Window is lower than 
that in observations. The predicted proper motions along 
the longitude in the three windows are greater than ob- 
served. There are also more proper motion data in small 
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Figure 15. Density contours from tlie input model (solid lines) 
and from orbits (dashed lines) in the x — y and y — z pianos for 
model 23. The left and right panes show the results by using 1000 
and 400 mass cells, respectively. 
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Figure 18. Velocity and velocity dispersion distribu- 
tion for model 23 (flp = 60 km kpc~^, disk mass 
Md = 1.0 X IO^^Mq and bar angle 6»bar = 20°). The solid 
line is for the model while the filled circles are data from 
BRAVA. 
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Figure 16. Projected density contours from the input model 
(solid lines) and from orbits (dashed lines) for model 23. 1000 
cells are used for solving the orbit weight. 
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fields from OGL E feumi et al ] |2004 ': 'Rattenburv et al."2007|; 
lRattenburv~M ao 2008,) and HST (Kozlowski et al. 200^, 
we do not compare with the proper motion in these field 
because the sky area in our model is 1° x 1°, much larger 
than the observed field size. 

4-3.2 Phase-space distribution and orbit families 

The phase-space distribution and the orbit family are use- 
ful to help us understand the model. Figure [19] shows the 
distribution of average energy versus angular momentum 
along the z-axis for non-zero weight orbits. The solid and 
dashed lines are the theoretical distributions of the energy 
versus angular momentum along the z-axis for retrograde 
and prograde motions from to 3kpc, respectively. The en- 
ergy and angular momentum along the 2-axis in the lab- 
oratory frame are defined as £;iab(r) = $(r) + [±|K|]V2, 
Jz,iab(7') — [±|T4|] X r. In the rotating frame, they are 
Erot{r) = <l>(r) + [±|K| + \np\rf/2, and J...ot(r) = [±|K| + 
Slp|r] x r. Here r is the radius, 14 is the circular velocity 
and $ is the potential of the system. The and '-' signs 
mean the retrograde/prograde motions. It is seen that most 
orbits are located in the range between the prograde and 
retrograde motions; only small number of orbits are disk 
orbits. For the orbit classifi cation, we use the method of 
ICarpintero fc Aguilaij (|l998l ). In Tabled we show the rel- 
ative fractions of orbit families with non-zero weights for 
model 23. It is seen that irregular orbits dominate the orbit 
families, the fraction of these orbits is over 90% for the mod- 
els if the integration time is one Hubble time. We also find 
that the relative orbit fraction strongly depends on the in- 
tegrated orbit time, the fraction of irregular orbits increases 
with increasing orbit integration time. 



Figure 17. Departure from self-consistency S as a function of the 
number of orbits. The star and diamond symbols represent the 
fitting results by using 1000 and 400 mass cells, respectively. 
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Figure 19. Distribution of the average energy versus angular 
momentum along the z-axis for orbits with non-zero weights. The 
upper and lower panels show the results in the laboratory and ro- 
tational frames, respectively. The filled circles, crosses and trian- 
gles represent orbits with large, intermediate and small weights, 
respectively. The solid and dashed lines represent the theoret- 
ical retrograde and prograde motions from to 3kpc, respec- 
tively. The bar model is model 23 with fip = 60 km s~^ kpc~^, 
Md = 5.0 X Mbar and Star = 20°. 



Table 2. for different input models, which are constrained by 
the velocity and velocity dispersion in four windows (fe = —4°, b = 
-6°,b = -8° and I = 0°). 
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Table 3. Observed proper motion dispersions in some fields. The 
bottom four rows are predictions from model 23. 



Field 


Q,bl 




(ji 




Ref. 




n 




(mas yr^-"-) 


(mas yr^-"-) 




Baade's Window 


(lr4) 




3.2 ±0.1 


2.8 ±0.1 


Spaenhauer et al. fl992) 


Baade's Window 


(1,-4) 




3.14 ±0.11 


2.74 ±0.08 


Zhao et al. (1996) 


Baade's Window 


(1.13,-3 


77) 


2.9 


2.5 


Kuiikcn & Rich (2002) 


Baade's Window 


(1,-4) 




2.87 ±0.08 


2.59 ± 0.08 


Kozlowski et al. (2006) 


Baade's Window 


(0.9,-4) 




3.06 ±0.11 


2.79 ± 0.13 


Soto et al. C2007~) 


Baade's Window 


(1,-4) 




3.13 ±0.16 


2.50 ±0.10 


Babusiaux et al. f2010) 


Baade's Window 


(1.13,-3 


76) 


3.11 ±0.08 


2.74 ±0.13 


Soto (2012) in preparation 


Plaut's Window 


(0,-8) 




3.39 ±0.11 


2.91 ±0.09 


Vieira et al. f2007. 2009) 


Sagittarius I 


(1.25,-2 


65) 


3.3 


2.7 


Kuiiken & Rich f2002) 


Sagittarius I 


(1.27,-2 


66) 


3.07 ±0.08 


2.73 ± 0.07 


Kozlowski et al. (2006) 


Sagittarius I 


(1.25,-2 


65) 


3.067 


2.760 


Clarkson et al. f2008) 


Sagittarius I 


(1.26,-2 


65) 


3.56 ±0.08 


2.87 ±0.08 


Soto (2012) in preparation 


NGC 6558 


(0.28,-6 


17) 


2.90 ±0.11 


2.87 ±0.13 


Soto (2012) in preparation 


Baade's Window 


(1,-4) 




4.44 


2.52 


Model 23 


Plaut's Window 


(0,-8) 




5.28 


2.32 


Model 23 


Sagittarius I 


(1,-3) 




4.43 


2.67 


Model 23 


NGC 6558 


(0,-6) 




4.46 


2.36 


Model 23 
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Table 4. Tangential velocity dispersion in some fields. The val- 
ues given in the references are derived from proper motions by 
assuming a distance to the Galactic centre Ro = 8kpc. 



Field 


(l,b) 










Rcf. 




( ) 




(km s 


) 


(km s ) 




Baade's Window 


(1,-4) 




121 ± 


4 


106 ±4 


Spacnhaucr ct aL (_1992) 


Baade's Window 


(1,-4) 




119 ± 


4 


104 ±3 


Zhao ct ah (1996) 


Baade's Window 


(1.13,-3 


77) 
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100 


Kuiiken & Rich (2002) 


Baade's Window 


(1,-4) 




109 it 


3 


98 i 3 


tvoziowsKi et al. izUUo) 


Baade's Window 


(0.9,-4) 




116 ± 


4 


106 ±5 


Soto et al. (2007') 


Baade's Window 


(1,-4) 




119 ± 


6 


95 ±4 


Babusiaux et al. f2010) 


Baade's Window 


(1.13,-3 


76) 


118 ± 


3 


104 ±5 


Soto et al. (2012) 


Plant's Window 


(0,-8) 




129 ± 


4 


110 ±4 


Vieira et al. (2007. 2009) 


Sagittarius I 


(1.25,-2 


65) 


123 




108 


Kuiiken & Rich (2002) 


Sagittarius I 


(1.27,-2 


66) 


117 ± 


3 


104 ±3 


Kozlowski ct al. (2006) 


Sagittarius I 


(1.25,-2 


65) 


116 




105 


Clarkson ct al. (2008) 


Sagittarius I 


(1.26,-2 


65) 


135 ± 


3 


109 ±3 


Soto et al. (2012) 


NGC 6558 


(0.28,-6 


17) 


110 ± 


4 


109 ± 5 


Soto et al. (2012) 


Baade's Window 


(1,-4) 




138 




90 


Model 23 


Plaut's Window 


(0,-8) 




121 




71 


Model 23 


Sagittarius I 


(1,-3) 




145 




97 


Model 23 


NGC 6558 


(0,-6) 




125 




78 


Model 23 


Baade's Window 


(1,-4) 




140 




106 


ZH96 


Sagittarius I 


(1,-3) 




146 




118 


ZH96 
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Table 5. Relative fractions of orbit families with non-zero weights 
for model 23. 



Orbit families 


Fractions" 


Fraction'' 


Fraction^ 


Fraction'* 


irregular 


92.9% 


91.2% 


81.3% 


80.6% 


open z tube 


5.8% 


5.6% 


13.5% 


12.0% 


thin z tube 


0.5% 


2.2% 


2.27% 


2.6% 


others'^ 


0.8% 


1.0% 


2.93% 


4.8% 



" The integration time of orbits is one Hubble time. 

^ The integration time of orbits is one half Hubble time. 

The integration time of orbits is one third Hubble time. 

The integration time of orbits is one fourth Hubble time. 
" others mean closed z tube, open box, open x tube etc. 
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5 STABILITY OF THE MODEL 

A perfect dynamical model must be self-consistent and sta- 
ble. From the projected density map and radial kinemat- 
ics, we know that our model is self-consistent. However, its 
dynamical stability still needs to be demonstrated. For a 
non-spherical model, N-body simulations can help us check 
the stability of the model . For this purp ose, we use the ver- 
satile code GADGET-2 (|Springell [2005h to run our galaxy 
simulation. 

One important step in N-body simulations is t o gen- 
erate the initial condition s. In our study, we follow IZhad 
l|l996t) and lWu et al] (|2009l '). The total number of initial par- 
ticles is 2 X 10^. They are randomly sampled from non-zero 
weight orbits. In each orbit, the selected number of particles 
is proportional to its orbit weight. The disk part is not in 
equilibrium outside 3 kpc, and 30% of the disk mass escaped 
during the simulation. The bar is essentially in equilibrium 
with an overall drift of 2 km s~^ perhaps due to the asym- 
metric escape of the disk particles. 

The N-body simulation for our best-fit model 23 is run 
for 2 Gyr, and 30 snapshots of particle positions and veloci- 
ties are stored. Figure [20l shows the snapshots of the model 
at 0, 0.1 and 0.5 Gyr. 

Two imperfections are noted: by construction, our disk 
part is not in equilibrium outside 3 kpc where no self- 
consistency was imposed in the Schwarzschild model. As a 
result about 30% of the disk mass evaporates due to the 
relatively shallow potential of the live particles. The bulk of 
the disk remains, however, the non-equilibrium of the disk 
makes it difficult for us to assess the long-term stability of 
the bar model. Secondly the center of the bar begins to drift 
slightly from the origin. The drift speed is very low, about 
2 km s~^, probably due to the recoil momentum of asym- 
metric evaporation of disk particles beyond 3 kpc. 

Nevertheless an examination of the moment of inertia 
IxY, Ixx and Iyy reveals that the bar rotates about 5 
times (i.e., there are 10 peaks in Ixv, Ixx and Iyy) in 0.5 
Gyrs, which is consistent with a constant pattern rotation 
speed 60 km s~^ kpc^^. The axis ratio appears stable as 
well with very little evolution within 2 Gyrs. Clearly more 
runs are necessary with particular attention to include a disk 
component. 



6 CONCLUSION AND DISCUSSION 

We have constructed 36 nuclear models using 
Schwarzschild's method, varying the bar angle, bar pattern 
speed and disk mass. Through fitting, we find the best-fit 
model has ftp = 60 kms"^kpc"\ Md = 1.0 x lO^M© and 
^bar ~ 20°. The model can reproduce the three-dimensional 
density, surface density and BRAVA velocity and velocity 
dispersion well. We tested two different methods of smooth- 
ing and two methods of generating initial conditions; our 
results are independent of these. 

Compared to the model of ZH96, our model can better 
reproduce the rotation curve as seen in the average radial 
velocity and the surface brightness distributions. However, 
our model is incomplete in at least two aspects, the first is 
that the predicted proper motions appear to be too high, 
and the second is the stability of the system is far from 
perfect. We discuss these two issues in turn. 
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Figure 20. Snapshots for the steady-state model of N-body sim- 
ulations for model 23. From left to right, the output time is 0, 0.1 
and 0.5 Gyr, respectively. The solid line at the bottom left panel 
indicates the line of sight to the Galactic centre. 

The proper motion dispersions in some windows pre- 
dicted from the model are higher than observed along the 
Galactic longitude. Possible reasons are: 

(i) The nuclear model in our paper is not perfect. Al- 
though we have calculated 36 models with different param- 
eters, it is possible that we still miss the model with the 
right combinations. In particular, we have used a fixed bar 
mass (2.0 x IO^'^Mq) in this study. Moreover, the density 
models adopted in this paper are obtained from fitting the 
surface brightness of the inner Galaxy. It has been shown 
that models with dif ferent axis ratio s can fit the surface 
brightness of the bar (|Zhao et al.lll99(j V In other words, the 
three-dimensional density distribution is not unique. Our 
model shows a stronger anisotropic distribution in proper 
motion than in observations, thus a less triaxial bar model 
may be preferred. 

(ii) The proper motion is obtained by using the orbits 
inner 5kpc around the Galactic center. If we use the orbits 
only inner 2.5 kpc around the Galactic center, the predicted 
value of proper motion will be changed. For example, in 
BW, ai = 3.72 mas yr~^ and CT(, — 2.53 mas yr~^, then the 
agreement between the model prediction and observation 
improves. A useful way to compare the model prediction 
with those observations is to use both the tangential velocity 
dispersion and proper motion dispersion. It is seen that the 
predicted velocity dispersion in the Plant's Window along 
the longitude direction is close to the observed value (See 
Table H]) . 

(iii) The predictions in our model are for pixels of 1° x 1°, 
the observed regions are far smaller. The predicted proper 
motions of model 23 in the range I = [—12°, 12°], b — 
[—10°, 10°] are available online q Future observations of 
large proper motion samples can help us to further constrain 
the models. 

The second shortcoming of our model is the dynami- 
cal instability. N-body simulations show that the model is 
only stable for 0.5 Gyr. The initial conditions for N-body 

^ http: / /cosmology.bao.ac.cn/" wangyg/ 
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Figure 21. Density contours from tlie input model (solid lines) 
and from orbits (dashed lines) in the x — y and y — z planes for 
model 23. The left panel is for the result of the orbit weights 
are solved by only using the three-dimensional density while the 
right one is that the orbit weights are solved by using the three- 
dimensional density, velocity and velocity dispersion. 

are generated from the orbit weights. In our model, most 
orbits are irregular. The fraction of irregular orbit strongly 
depends on the potential. Our model includes a prolate bar, 
a boxy bulge and an axi-symmetric disk; the axis ratios of 
these three components are different. In particular, the pres- 
ence of the bar implies the axi-symmetric disk should not be 
present in the central part since no circular orbits can exist. 
This may have limited the dynamical stability of our sys- 
tem. In addition, we only consider the self-consistency of the 
model within 3 kpc. The disk is an important component in 
our model which dominates the mass beyond 3 kpc. We also 
check the self-consistency of the model inside 6 kpc, which 
covers the outer Lindblad resonance region (^ 6kpc). Fig- 
ure [21] compares the input and the reconstructed densities 
from orbits. As can be seen, the agreement is good within 
3 kpc, but the scatters become increasingly larger (~ 30%) 
beyond 3 kpc. The N-body simulation again shows that the 
model is only stable within 0.5 Gyr. In the future, it may be 
desirable to start with a bar model f rom N-body simulations 
such as that from IShen et all (|2010l ). 

At present, most density or potential models of the 
Galactic nucleus are constructed using photometric data 
(surface brightness and star counts) alone. However, no den- 
sity model is constructed including the kinematics, such as 
velocity, velocity dispersion and proper motion. In the fu- 
ture, it may be useful to construct the Galactic bar model 
using the density and photometry at the same time. Future 
observations such as GAIA will be particularly valuable to 
help us construct a density model. 
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